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Abstract 


The  nonlinear  composite  analysis  code  “LAMP AT”  is  implemented  in  the  nonlinear  explicit 
finite  element  code  LLNL-DYNA3D  as  a  new  user-defined  material  model.  All  subroutines  of 
LAMP  AT  are  implemented  as  user-defined  Material  Model  46.  The  user-defined  material 
subroutine  calls  an  external  data  base  file  that  contains  material  properties  for  several 
composites.  In  addition,  LAMP  AT  is  modified  for  use  in  an  explicit  time-integration  solver. 

The  model  is  improved  to  account  for  loss  of  symmetry  of  the  material  stiffness  matrix  resulting 
from  degradation  of  the  elastic  moduli  during  damage  evolution.  The  model  implementation  is 
validated  through  a  one-element  simulation  and  penetration  simulations.  In  addition,  comparing 
the  prediction  of  the  LSDYNA  elastic  material  model  with  that  of  LAMP  AT  validates  the 
implementation. 
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1.  Introduction 


The  U.S.  Army  Research  Laboratory  (ARL)  is  actively  seeking  better  computational  material 
models  to  assist  the  design  and  optimization  process  of  armors.  The  composite  lightweight 
integral  armor  under  consideration  contains  a  variety  of  materials,  such  as  fiber-reinforced 
composites,  metals,  loss  fabric,  ceramics,  and  rubber.  To  be  able  to  simulate  penetration  and 
perforation,  a  robust  computer  code  must  be  utilized.  The  nonlinear  explicit  finite  element  code 
LLNL-DYNA3D*  [1]  is  an  excellent  candidate  to  successfully  simulate  such  problems. 

However,  the  code  is  very  general  and  needs  to  be  augmented  with  user-defined  material  models. 
The  user-defined  material  model  options  in  DYNA3D  allow  users  to  customize  the  code  with 
advanced  material  models  of  interest  to  the  armor  community.  Researchers  at  ARL  have  been 
developing  the  “LAMP AT”  nonlinear  composite  analysis  code  for  the  analysis  of  thick-section 
composite  structures  for  more  than  10  years  [2].  LAMP  AT  models  the  material  nonlinear 
behavior  observed  in  some  composites  and  is  based  on  homogenization  procedures.  To  be  able 
to  design  armored  vehicles,  one  must  be  able  to  simulate  the  high-speed  ballistic  impact  event 
with  good  accuracy.  In  this  manner,  one  can  use  numerical  simulation  to  optimize  an  armor  that 
is  able  to  contain  high-speed  projectiles.  The  material  model  developed  in-house  at  ARL  can 
perform  this  task  provided  that  it  is  implemented  in  a  wave  propagation  code  like  DYNA3D. 

The  LAMP  AT  code  was  implemented  in  an  initial  effort  [3]  into  the  commercial  code  LSDYNA 
[4]  as  a  user-defined  material  subroutine.  However,  the  source  code  for  the  commercial  code  is 
not  available  for  further  development,  which  makes  it  more  difficult  to  link  to  more  sophisticated 
armor  optimization  algorithms  [5], 

The  implementation  of  the  LAMP  AT  code  into  the  DYNA3D  code  herein  is  described.  The 
modification  of  the  nonlinear  LAMP  AT  code  is  necessary  for  successful  utilization  of  the  code 
in  the  finite  element  method.  This  modification  is  described  next.  The  implementation  is 
validated  through  several  examples  that  are  described  and  presented  in  the  following  sections. 

A  listing  of  the  modified  LAMP  AT  code  is  presented  in  Appendix  A.  To  understand  the  solid 
element  formulation,  Appendix  B  lists  the  algorithm  that  completely  describes  this  element. 


2.  Implementation 


The  original  LAMP  AT  development  is  performed  for  analysis  of  composite  structures  with  a 
mechanics  of  composite  materials  approach.  In  that  approach,  the  stresses  are  incrementally 


LLNL-DYNA3D  is  the  Lawrence  Livermore  National  Laboratories  version  of  the  finite  element  computer  code  hereafter 
referred  to  as  simply  DYNA3D. 
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increased,  and  the  corresponding  strains  are  obtained  (or  vice  versa).  The  original  formulation 
leads  to  numerical  instability  due  to  loss  of  symmetry  of  the  material  stiffness  matrix  when 
implemented  into  three-dimensional  finite  element  codes  like  DYNA3D.  The  loss  of  symmetry 
in  the  material  stiffness  matrix  is  attributed  to  how  the  elastic  constants  are  degraded.  The 
material  reciprocity  relation  for  an  orthotropic  material  is  given  by 
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The  compliance  matrix  can  lose  symmetry  if  material  constants  are  arbitrarily  degraded  during 
damage  progression.  The  loss  of  symmetry  is  avoided  by  utilizing  the  following  compliance 
matrix: 
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The  material  nonlinearity  treatment  in  explicit  finite  elements  is  different  than  that  for  implicit 
finite  elements.  In  implicit  finite  elements,  we  have  load  increments  and  equilibrium  iterations. 
The  tangential  stiffness,  due  to  material  nonlinearity,  must  be  updated  at  each  load  increment  and 
should  also  be  updated  at  each  equilibrium  iteration.  The  material  nonlinear  behavior  in 
LAMP  AT  is  summarized  by  the  following  equation: 
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The  concept  of  equilibrium  is  given  by  equation  5.  When  the  internal  force  vector  is  equal  to  the 
external  force  vector,  equilibrium  is  satisfied.  The  internal  force  vector  is  the  integral  over  the 
volume  of  the  stresses  times  the  strain  displacement  matrix  of  the  finite  element.  The  internal 
force  vector  is  defined  by  equation  6. 


Rml  =/?“'. 

R*  =  \BTCydv. 

In  implicit  finite  elements,  the  external  load  is  incremented  in  several  steps.  At  the  first  iteration 
of  the  first  step,  the  internal  forces  are  normally  not  in  equilibrium  with  the  external  forces.  The 
out-of-balance  force  is  called  the  residual  A R .  Due  to  this  residual,  a  correction  to  the 
displacement  is  obtained  by  inverting  K  in  equation  7  and  solving  for  the  displacement. 


(5) 

(6) 


RM-Ra'  =AR  =  [ZT]{Au}.  (7) 

Aw  =>  Ae  =>  Act  .  (8) 

From  the  corrected  displacement,  a  new  strain  increment  is  obtained,  and  a  new  stress  increment 
can  be  calculated  (equation  8).  Subsequently,  new  internal  forces  are  obtained.  These  iterations 
are  continued  until  the  out-of-balance  force  is  less  than  a  user-defined  value. 

The  iteration  algorithm  previously  described  can  be  understood  by  examining  Figures  1  and  2. 

At  the  start  of  the  load  step  “n”  there  is  a  tangential  material  stiffness  update.  As  equilibrium 
iterations  proceed,  the  material  stiffness  matrix  is  also  updated  as  depicted  in  Figure  2. 

In  explicit  finite  elements,  however,  the  load  varies  incrementally  at  each  time  step.  There  is  no 
equilibrium  iteration,  since  the  acceleration  is  solved  for  in  the  same  time  step.  The  acceleration 
is  integrated  once  to  obtain  velocity  and  once  again  to  obtain  displacement.  In  explicit  finite 
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elements,  we  have  many  load  increments  due  to  many  cycles  (time  steps)  to  maintain  time 
integration  stability  that  is  governed  by  the  Courant  condition.  The  integration  time  steps  can  be 
on  the  order  of  microseconds  or  fractions  of  microseconds.  Consequently,  in  a  typical  impact 
simulation  thousands  of  cycles  are  possible.  This  fact  can  render  the  frequent  update  of  the 
tangential  material  stiffiiess  matrix  computationally  inefficient.  LAMP  AT  calculations  are 
computationally  costly  due  to  several  transformations  of  stresses  and  strains  for  multiple  layers 
from  global  to  local  and  back  from  local  to  global  frames  of  reference.  In  explicit  finite 
elements,  the  LAMP  AT  material  stiffiiess  matrix  updates  must  be  limited  to  achieve 
computational  efficiency  for  large-scale  finite  element  impact  simulations.  For  this  reason,  the 
tangential  material  stiffiiess  matrix  is  stored  in  history  variables  for  subsequent  use  in  the 
calculation.  The  current  implementation  permits  the  user  to  define  the  frequency  of  the  material 
stiffiiess  matrix  update. 

The  initial  verification  of  the  implementation  is  performed  on  an  isotropic  elastic  material.  A 
linear  elastic  material  behavior  is  obtained  in  LAMP  AT  by  suppressing  the  power  law  stress- 
strain  relation  (equation  4).  The  material  model  is  then  tested  for  prediction  of  stress  for  an 
isotropic  elastic  material  by  keeping  the  composite  material  transformation  matrices  intact. 

Exact  prediction  is  obtained  as  compared  to  the  existing  material  elastic  material  model  in 
DYNA3D.  Finally,  the  model  is  tested  for  material  nonlinear  behavior. 


3.  Results 


A  finite  element  model  of  a  block  consisting  of  one  element  is  developed  for  verification  of  the 
nonlinear  behavior.  The  block  is  constrained  at  one  end  and  a  prescribed  displacement  is  applied 
to  the  other  end.  Two  material  systems  are  considered  for  validation  of  the  implementation.  An 
isotropic  elastic  material  (aluminum)  and  nonlinear  elastic  composite  (graphite/epoxy)  are 
considered.  The  linear  elastic  material  is  not  shown  here,  as  it  leads  to  the  same  behavior  of  the 
original  elastic  material  model  in  DYNA3D.  Figures  3-5  depict  the  stress  vs.  strain  for  the 
material  system  considered.  The  validation  of  the  nonlinear  behavior  is  carried  out  in  the  three 
coordinate  directions  (i.e.,  x-,  y-,  and  z-directions).  The  nonlinear  material  behavior  is  illustrated 
in  these  plots. 

Next,  the  implementation  is  tested  for  numerical  stability  and  the  ability  to  simulate  an  impact 
event.  For  this  purpose,  two  models  are  considered.  The  first  model  is  a  ball  impacting  thick 
plates  (Figures  6  and  7).  The  plates  are  made  of  graphite/epoxy  with  nonlinear  material 
behavior.  The  impact  velocity  is  1.5e+4  in/s.  The  second  example  considered  is  a  flat  plate 
impacting  a  curved  plate.  The  flat  plate  is  given  two  velocity  components.  One  component 
(2.0e+4  in/s)  is  in  the  normal  direction  and  the  other  velocity  component  (1.0e+4  in/s)  is  in  the 
circumferential  direction.  The  model  and  the  result  of  the  simulation  are  depicted  in  Figures  8 
and  9.  The  flat  plate  is  made  of  graphite/epoxy  with  nonlinear  material  behavior. 
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Finally,  a  12  x  12  *  1.6-in  composite  plate  made  of  S-2  Glass  /epoxy  was  constructed,  clamped 
at  four  comers,  and  impacted  with  a  .50-cal.  (1.6-oz)  trimetallic  bullet  consisting  of  a  copper 
jacket,  steel  core,  and  lead-tipped  filler.  The  experimentally  determined  V50  velocity  was 
17,280  in/s  (1440  fps),  although  this  value  was  not  known  when  the  simulations  were  performed. 
Three  different  composite  material  models  were  evaluated  to  see  which  could  best  predict  the 
V50  that  was  determined  experimentally.  The  first  model  is  the  orthotropic  elastic  Material 
Model  22  in  LSDYNA.  The  second  model  is  the  micromechanical  strain-rate-sensitive 
unidirectional  composite  material  model  described  by  Tabiei  and  Chen  [6]  and  implemented  as  a 
user-defined  material  subroutine  in  LSDYNA.  The  third  material  model  is  the  new 
LAMP  AT/D  YNA3D  implementation.  The  prediction  of  each  of  the  material  models  is 
considered  and  described  in  the  next  section. 

1.1  Material  22  in  LSDYNA 

The  composite  plate  is  modeled  with  45  layers  of  solid  elements  through  the  thickness.  Each 
layer  represents  a  layer  of  unidirectional  composite.  Two  layers  together,  one  in  the  zero 
direction  and  the  other  in  the  90  direction,  represent  one  layer  of  plain  weave  composite.  The 
material  (S-2  Glass/epoxy)  properties  for  Material  Type  22  used  in  the  simulation  are  listed  in 
Table  1.  Several  simulations  were  carried  out  with  bullet  velocities  that  ranged  from  13,000  to 
16,000  in/s.  Figures  10-15  depict  the  results  for  varying  the  impact  velocity  to  find  the  V50. 
Figure  10  shows  the  penetration  of  the  bullet  into  the  plate,  the  velocity  of  the  bullet  as  a 
function  of  time,  and  the  displacement  of  the  center  of  the  plate  as  a  function  of  time.  The 
velocity  of  node  52569  is  shown  in  Figure  10.  This  node  is  on  the  center  of  the  back  of  the 
bullet.  The  displacement  of  node  1  is  also  presented  in  Figure  10.  This  node  is  on  the  center 
back  face  of  the  plate.  It  can  be  seen  that  the  bullet  is  stopped  in  the  plate  with  an  impact 
velocity  of  13,000  in/s.  Figures  1 1  and  12  depict  results  for  impact  velocities  of  13,500  and 
14,000  in/s,  respectively.  In  both  of  these  simulations,  the  bullet  stopped  at  about  250  ms.  The 
maximum  deflection  of  node  1  in  both  cases  is  around  0.03  in.  Note  that  for  the  case  of  the 
14,000  in/s  impact  velocity,  the  maximum  deflection  is  0.14  in.  This  is  due  to  the  eroding 
algorithm.  When  an  element  fails,  it  is  removed  from  the  database.  The  nodes  of  the  removed 
element  are  also  removed.  The  sharp  increase  of  the  displacement  in  Figure  12  has  no  physical 
meaning.  Figure  13  depicts  the  V50  velocity  results  (14,500  in/s).  The  impact  velocity  is 
increased  to  15,000  and  16,000.  The  results  for  these  cases  are  shown  in  Figures  14  and  15.  It 
can  be  concluded  that  velocities  greater  than  14,500  in/s  lead  to  perforation  of  the  bullet  with  a 
non-zero  residual  velocity.  The  predicted  V50  velocity  is  between  14,000  and  14,500  and 
underpredicts  the  experimental  value  by  16%-19%.  The  dynamic  deflection  of  the  center  of  the 
plate  for  this  case  is  ~0.05  in. 


*  S-2  Glass  is  a  registered  trademark  of  Owens  Coming  Corp. 
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1.2  Micromechanical  Material  Model  in  LSDYNA 

Tabiei  and  Chen  developed  a  micromechanical  material  model  for  unidirectional  composites  [6]. 
The  material  model  assumes  linear  elastic  behavior  up  to  the  first  constituent  failure. 
Microfailure  criteria  are  utilized  in  the  model.  Tabiei  et.  al  [7]  also  presented  a  nonlinear 
viscoplastic  strain-rate-sensitive  micromechanical  material  model  with  progressive  microfailure 
criteria.  Both  models  are  programmed  as  user-defined  material  routines  in  the  explicit  code 
LSDYNA.  The  material  models  are  used  to  simulate  the  bullet  penetration  problem.  The 
material  (S-2  Glass/epoxy)  properties  used  in  the  simulation  are  listed  in  Table  2.  Note  that  in 
this  material  model,  the  properties  of  each  of  the  constituents  (fiber  and  resin)  must  be  provided. 
Figure  16  shows  the  prediction  of  the  V50  velocity  of  about  14,250  in/s,  which  underpredicts  the 
experimental  value  by  ~17.5%.  Several  computational  iterations  were  performed  to  obtain  this 
velocity,  but  are  not  illustrated  in  this  report.  Two  sets  of  simulations  are  shown  in  Figure  16; 
Figure  16a  is  the  simulations  obtained  with  the  strain-rate-sensitive  micromechanical  model,  and 
Figure  16b  is  the  simulations  obtained  with  the  strain-rate-insensitive  micromechanical  model. 
The  predictions  of  the  rate  insensitive  model  are  stiffer  than  those  for  the  rate  sensitive  model; 
this  can  be  seen  from  the  resultant  velocity,  the  displacement  of  the  bullet,  and  the  plate  dynamic 
deflection.  The  increased  stiffness  observed  for  the  rate  insensitive  micromechanical  model  is 
due  to  the  nonlinear  material  behavior  and  strain  softening  for  the  matrix  material.  The  dynamic 
deflection  of  the  plate  is  predicted  to  be  about  0.15  in,  which  is  larger  than  the  LSDYNA 
Material  Model  22  prediction  of  0.05  in.  In  both  models,  the  effective  strain  is  used  as  an 
ultimate  criterion  for  element  erosion.  The  value  of  the  effective  strain  at  failure  is  taken  to  be 
20%  for  both  cases.  However,  in  the  future,  a  strain-rate-dependent  failure  criterion  might  be 
employed.  Note  that  the  prediction  of  V50  for  the  rate  sensitive  material  model  is  less  than  the 
prediction  of  the  V50  with  no  strain  rate  effect.  This  attributed  to  the  fact  that  the  strain-rate- 
sensitive  material  model  considers  material  nonlinear  behavior  in  the  composites  (strain 
softening).  However,  the  Material  Model  with  no  rate  effect  assumes  that  everything  is  elastic 
up  to  the  point  of  failure.  The  principal  author  will  be  implementing  a  Cooper-Symond  type  of 
rate  sensitive  failure  criterion  in  the  referenced  micromechanical  model. 

13  LAMPAT/DYNA3D  Implementation 

After  validating  the  formulation  for  nonlinear  behavior  (Figures  3-5)  the  implemented  material 
model  is  utilized  to  simulate  the  bullet  penetration  problem.  The  aim  is  to  predict  the  V50 
velocity  using  the  LAMPAT/DYNA3D  code.  The  input  file  for  the  bullet  model  concerning  the 
data  required  for  LAMP  AT  material  subroutine  is  listed  in  Table  3.  The  material  type  is 
programmed  as  Material  46.  The  third  row/first  column  requires  an  input  of  the  elastic  modulus. 
The  fourth  row/first  column  requires  the  Poisson’s  ratio.  These  two  material  constants  are  used 
in  the  time  step  calculation  and  the  penalty  contact  stiffiiess  calculation  by  DYNA3D.  The  fifth 
row/first  column  requires  the  interval  frequency  for  material  stiffiiess  updates.  Finally,  the  sixth 
row/first  column  requires  the  allowable  for  effective  strain  at  failure  for  the  eroding  algorithm  in 
DYNA3D.  The  same  value  (20%  failure  strain)  is  used  for  eroding  as  the  other  two  previously 
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presented  material  models.  Table  4  lists  the  material  properties  used  by  LAMP  AT/D  YNA3D  for 
the  S-2  Glass/epoxy.  There  is  difficulty  defining  the  exponential  parameters  for  the  material 
nonlinear  behavior  in  LAMP  AT  since  they  have  not  been  determined  experimentally,  but  for  the 
purpose  of  comparison,  the  same  exponents  used  in  the  graphite/epoxy  material  [2]  are  used  for 
the  S-2  Glass/epoxy.  Two  extreme  cases  are  considered  for  the  simulation  of  the  bullet 
penetration  problem.  In  the  first  simulation  (case  1)  the  composite  plate  is  assumed  to  be  linear 
up-to-failure  in  LAMP  AT;  this  behavior  is  achieved  by  suppressing  the  nonlinear  update  of  the 
elastic  and  shear  modulus  in  LAMP  AT.  In  the  second  simulation  (case  2),  nonlinear  material 
behavior  is  assumed  in  LAMP  AT.  Figures  17  and  18  show  the  penetration  predictions  of  the 
two  cases  considered.  In  both  cases,  the  material  stiffness  matrix  is  updated  every  100  cycles. 
Figure  17  is  for  linear  up-to-failure  material  behavior,  and  Figure  18  is  for  nonlinear  up-to- 
failure  material  behavior.  It  is  apparent  that  there  is  significant  difference  in  the  results.  The 
nonlinear  material  behavior  (case  2)  predicts  a  much  softer  response  than  the  prediction  of  the 
linear  material  behavior  (case  1).  For  case  2,  the  bullet  is  almost  intact,  and  for  case  1  the  bullet 
head  and  jacket  is  fully  eroded.  A  comparison  is  also  made  for  the  effect  of  the  interval  of 
material  stiffness  matrix  updates  on  the  behavior  of  the  considered  problem.  Figures  19  and  20 
show  the  effective  plastic  strain  in  the  bullet  and  the  effective  stress  in  the  plate  for  two 
frequencies  of  the  material  stiffness  matrix  updates. 

Figure  19  shows  the  result  for  every  10  cycles  update,  and  Figure  20  shows  the  result  for  every 
100  cycles  update.  For  the  10  cycles  update,  the  effective  plastic  strain  in  the  bullet  is  predicted 
to  be  0.744,  and  the  effective  stress  in  the  plate  is  245  kpsi  at  60  ms.  For  the  100  cycles  update, 
the  effective  plastic  strain  in  the  bullet  is  predicted  to  be  0.791,  and  the  effective  stress  in  the 
plate  is  457  kpsi  at  60  ms.  The  CPU  time  difference  between  the  two  cases  is  almost  1  order  of 
magnitude.  This  makes  the  analysis  with  few  cycles  for  material  stiffness  update  uneconomical 
even  for  small-scale  simulations. 


4.  Summary 


The  LAMP  AT  nonlinear  composite  material  code  is  successfully  implemented  in  the  nonlinear 
dynamic  explicit  finite  element  code  LLNL-DYNA3D.  LAMP  AT  is  added  as  a  new  material 
model  in  DYNA3D  and  is  modified  to  fit  suitably  in  an  explicit  time  integration  scheme.  The 
code  calls  an  external  database  file  that  contains  properties  of  several  composite  materials. 
Therefore,  most  of  the  material  data  is  read  from  the  external  file,  and  the  limitation  on  the 
number  of  material  constant  by  DYNA3D  is  eliminated.  The  implementation  is  verified  for 
linear  and  nonlinear  behavior  of  composites.  The  nonlinear  behavior  is  verified  by  considering  a 
block  of  material  that  is  loaded  in  all  three  coordinate  directions.  Impact  simulations  are  also 
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carried  out  to  investigate  the  stability  of  the  LAMP  AT  code  in  the  explicit  finite  element  code 
DYNA3D.  A  bullet  penetration  problem  is  considered  and  several  simulations  are  carried  out 
using  three  material  models.  Material  type  22  (laminated  composite)  in  LSDYNA,  a  strain-rate- 
sensitive  micromechanical  model,  and  LAMPAT/DYNA3D  are  considered,  and  their  predictions 
are  presented.  The  following  can  be  summarized  from  the  presented  simulation  of  the  bullet 
penetration  problem: 

•  Material  Model  22  in  LSDYNA  underpredicted  the  V50  velocity  by  16%-19%;  however, 
the  prediction  of  dynamic  deflection  is  very  small  (0.03  in). 

•  The  micromechanical  material  model  predicted  the  V50  velocity  comparable  to  Material 
Model  22  in  LSDYNA.  The  dynamic  deflection  prediction  is  significantly  higher  than 
Material  Model  22  (0.15  in)  and  seems  more  reasonable,  although  no  experimental  data  on 
dynamic  deflection  were  experimentally  obtained. 

•  Strain-rate-sensitive  and  insensitive  micromodels  are  considered  for  the  bullet  penetration 
problem.  The  strain-rate-sensitive  material  model  considers  material  softening  and 
therefore  predicted  a  softer  behavior  than  the  rate  insensitive  model. 

•  The  LAMP  AT  arterial  nonlinear  input  data  cannot  easily  be  obtained;  therefore,  this 
presents  difficulty  for  simulating  penetration  experiments.  If  material  exponents  similar  to 
those  in  the  graphite/epoxy  are  used  in  the  S-2  Glass/epoxy,  premature  failure  in  the 
element  occurs.  This  leads  to  significantly  underpredicting  the  V50  velocity. 

•  The  CPU  time  required  by  Material  Type  22  in  LSDYNA  and  the  micromechanics  material 
model  is  significantly  less  than  that  of  LAMP  AT.  The  CPU  time  is  on  the  order  of  hours 
for  the  Material  22  and  the  micromechanics  material  model.  However,  it  is  on  the  order  of 
days  for  the  LAMPAT/DYNA3D  code. 

•  The  CPU  time  requirement  for  LAMPAT/DYNA3D  increases  1  order  of  magnitude  when 
the  frequency  of  the  material  stiffness  matrix  updates  is  decreased  by  1  order  of  magnitude. 

•  There  is  a  significant  difference  in  prediction  between  the  low-  and  high-frequency 
material  stiffness  matrix  updates.  Therefore,  for  more  accurate  predictions,  the  material 
stiffness  matrix  must  be  updated  more  frequently.  This  will  render  analysis  of  even  small 
finite  element  models  computationally  inefficient. 

The  LAMPAT/DYNA3D  code  could  be  improved  for  more  accuracy  in  predicting  impact  and 
penetration  problems.  Several  suggestions  are  presented  to  make  the  code  more  suitable  for 
finite  element  impact  simulations. 
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5.  Possible  LAMPAT  Improvements 


The  following  is  a  list  of  possible  improvements  that  could  make  the  LAMP  AT/D YNA3D 
material  model  more  accurate  and  robust  for  finite  element  impact  simulation  problems: 

•  Adding  a  prefailure  unloading  scheme.  The  current  implementation  assumes  that 
unloading  is  along  the  loading  path. 

•  Adding  a  postfailure  unloading  scheme.  The  current  implementation  also  assumes  that 
unloading  is  along  the  initial  loading  path. 

•  Adding  strain  rate  sensitivity  to  the  model. 

•  Adding  a  more  realistic  progressive  failure  process  in  the  homogenization  process. 

•  Postprocessing  can  be  improved  by  identifying  the  various  progressive  failure  or  damage 
mechanisms  as  they  evolve  temporally. 
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Figure  7.  Penetration  of  a  rigid  ball  into  three  composite  plates  (final  geometry). 


Figure  8.  Penetration  of  a  rigid  plate  into  a  composite  curved  plate  (initial  geometry). 
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Figure  9.  Penetration  of  a  rigid  plate  into  a  composite  curved  plate 
(intermediate  geometry). 
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Figure  14.  Impact  velocity  is  15,000  in/s. 
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Figure  15.  Impact  velocity  is  16. 
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Figure  17.  Linear  material  behavior  up-to-failure. 
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Figure  18.  Nonlinear  material  behavior  up-to-failure. 
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Figure  19.  Stiffness  update  every  10  cycles. 
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Figure  20.  Stiffiiess  update  every  100  cycles. 
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Table  1.  Material  22  in  the  LSDYNA  code. 
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Table  2.  The  strain-rate-sensitive  micromodel  as  user-defined  routine  in  LSDYNA. 
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Table  3.  The  LAMP  AT  material  model  in  DYNA3D. 


1  46  1.66800-4  0  0  0.0000000  0  0.0000000  0.0000000  000 
PART  PID  =  1,  COMPOSITE  0 

3.00000+7  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000 

0.3000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000 

1.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000 

1 . 00000e2  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000 

0.2000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000 

0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000  0.0000000 


Table  4.  The  data  file  for  the  S-2  Glass  in  LAMP  AT. 
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Appendix  A.  Brick  Element  Formulation  in  DYNA3D 


This  appendix  appears  in  its  original  form,  without  editorial  change. 
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Brick  Element  Formulation  in  DYNA3D 

The  following  is  a  brief  description  of  the  explicit  time  integration  algorithm  in  DYNA3D  for 
Solid  Element.  This  algorithm  is  reported  in  here  to  indicate  the  hierarchy  and  position  of 
material  model  in  the  solid  element  formulation. 


1 .  Set  initial  conditions  at  time  t=0: 

n  =  0;  At  =  0(timestep);X°  =  X(t  =  0);u°  =  u(t  =  0);  v"1/2  =  v(t  =  0). 
where  n  is  the  number  of  iteration  (number  of  cycle). 

2.  Go  to  step  3.5. 

3.  Explicit  time  integration 

3.1  Calculate  nodal  accelerations 

«-  =-^fcr' +F*' ~F* 

JW|Y 

where  is  the  nodal  mass  of  node  i.  And  C  is  the  damping  matrix. 

3.2  Impose  boundary  conditions 

3.3  Impose  nodal  constraints 

3.4  Update  velocities  and  geometry 

v-,/2  =v"-3/2  +a”-'At 

u"  =w"-1  +v"'1/2A/ 


Xn=X°+un 


3.5  Compute  external  nodal  forces  F£, 

■  evaluate  the  load  curves 

■  compute  surface  tractions 

■  compute  concentrated  loads 

■  compute  body  force  loads 

3.6  Process  solid  elements 

3.6.1  Compute  velocity  strains  and  spin  rates 
a).  Compute  Hourglass  Type  1, 2  or  4 

■  Compute  the  Jacobian  matrix  J(0, 0, 0) ,  the  element  volume  Ve,  and 
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00* 

dx 

H 

=  J-\ 0,0,0) 

dh 

dy 
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d*k 
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_  dz  _ 
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(k=l,2, 3,4) 
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■  Compute  the  velocity  gradient  matrix  L,  velocity  strains  s  and  spin  rates  co: 


/5v"_1/2  8 
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.  dr .  , 

V  *  /(o,o,o) 


v^y 


(0,0,0) 


+  [(V,4)-I,! -(V,6)— '=]! 


'V 


(0,0,0) 


y  (o.o.o) 

\  y  (o.o.o) 


#  #i-l/2  , 

f  =i[i”-i/2 +(r-,/2>r] 


(i,y=  1,2,3) 


»-l/2  1 


GT  ‘  ‘  =^[I"-1/2  -(X"“,/2)r] 


fl)-,,2A/  =  y[I-,/2  -(I"-,/2)r] 

■  Compute  efi  for  the  timestep  calculation. 

b) .  Compute  Hourglass  Type  3  or  5 

3.6.2  Compute  the  characteristic  length  le  for  the  timestep  calculation 

•--T- 

where  Vc  is  the  element  volume,  Atma  is  the  area  of  the  largest  side. 

3.6.3  Constitutive  Modeling 

a) .  Get  the  global  stresses  a^,  from  storage 

b) .  Compute  [a”-1  -  crn_1  (con"1/2  At)  +  (con~1/2 At)CTn~' 

c) .  Stress  calculation 
■  Find  local  coordinate  system  and  transformation  matrix. 

Find  the  unit  vectors  of  local  coordinate  system: 

_  r31  X  r42  .  _  r21  ~(r21  ‘e3)e3  .  _ 

|r31  X  r42  |  ||r21  ~  (r21  *  e3  )e3 1| 

where  rjj  denotes  the  vector  from  Node  j  to  Node  i. 

Assume  ei  =  l;i  +  mj  +  njk ,  and  i,  j  ,k  are  the  unit  vectors  of  global  system, 


J  global 
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we  have  the  transformation  matrix  T: 

/j  mx  «, 

L  tm,  «, 


r-[e,  8,  ej  = 


/3  m3  n3 


1  Transform  the  global  stresses  [or” 1  -  a"  '(to”  1/2At)  +  (o”  1/2At)crn  1  ]gjobal  to  local 
system: 

[<T’-'  ~<T-'  A<)ff  ]„ 


.  n-1/2 


1  Transform  the  global  velocity  strains  s  to  local  system: 


n-V2\ 

£ 

V  ) 


=  71 


local 


(  mn- 1/2  > 
£ 

V  y 


global 


■  Compute  the  value  of  max {C, , ,  C22 ,  C33 }  from  the  material  matrix  for  the  timestep 

calculation. 

■  Update  stresses  in  local  system 


k'L  =K~‘  (<'^)+(<,'^)<lw+cSB 


A t 


local 


■  Transform  the  stresses  cf  from  local  system  to  global  system 

^global  =  ^  local  T) 

d) .  Timestep  calculation 

■  Compute  the  adiabatic  sound  speed 

c  =  ^|max^n»^22^33}  ^  where  p  is  the  material  density. 

■  Compute  q  which  is  a  function  of  the  bulk  viscosity  coefficients  Co  and  Ci : 

f  C,c  +  C0/e|4|  for  eu  <  0 
9  [  0  for  4  >0 

where  /e  is  the  characteristic  length  of  the  solid  element. 

■  Compute  the  timestep  for  each  element 

le 

~  _  ,  /  2  .  2  \  1/2 
q  +  (q  +  c  ) 

and  take  the  minimum  value 

At  =  min {A^ ,  At2 AtN  } ,  where  N  is  the  number  of  elements. 

e) .  Energy  calculation 

f) .  Put  the  global  stresses  c”loba]  in  storage 

3.6.4  Calculation  of  Hourglass  forces  and  internal  forces  (F^  -F"t) 

■  Calculate  Hourglass  forces  for  elements  fh” 

■  Calculate  nodal  internal  forces  of  elements  f". 
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f£t  =  (BTadV  =  BT(0,0,0)anVe 

■  Assemble  (fhng  -  )  into  global  nodal  force  vector  (F£t  +  Fh"g  -  F", ) 

3.7  Prepare  to  start  the  next  step 

n  =  n+  1 

3.8  Go  to  step  3.1 

The  above  is  a  detailed  flow  chart  of  all  steps  necessary  to  understand  and  program  the 
LAMP  AT  within  the  solid  element.  As  can  be  seen  from  the  above  section  (or  bullet)  3.6.3  is 
where  the  modification  for  the  material  model  enters  the  formulation. 
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Appendix  B.  LAMPAT  FORTRAN  Source  Code 


This  appendix  appears  in  its  original  form,  without  editorial  change. 
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Intentionally  left  blank. 
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LAMPAT/DYNA3D  FORTRAN  CODE 


subroutine  f3dm46  (cm,bqs,nhex) 
parameter (lnv=32) 

implicit  double  precision  (a-h,o-z) 

LAMPAT  IN  DYNA3D,  By:  A.  TABIEI  (ala.tabiei@uc.edu),  Feb.  2001 
common/bkO  2 /dt 1 , dt2 , iburn, isdo, iorder 

common/bk06/time (2 , 8) ,head(12) , idmmy, iadd, ifil,maxsiz,ncycle 
common/bk20/idummy (10) ,ndum 

common/ aux2/dl  (lnv)  , d2  (lnv)  ,  d3  (lnv)  ,  d4  (lnv)  ,  d5  (lnv)  ,  d6  (lnv)  , 

1  wzzdt  (lnv)  ,  wyydt  (lnv)  ,  wxxdt  (lnv)  ,  einc  (lnv) 
common/ auxl4 / 

+  sigl (lnv)  ,  sig2 (lnv) , sig3 (lnv) , sig4 (lnv) , 

+  sig5 (lnv) , sig6 (lnv) , epx (lnv) , davg (lnv) , 

+  hisvl (lnv) , hisv2 (lnv) , hisv3 (lnv) ,hisv4 (lnv) ,hisv5 (lnv) , 

+  hisv6(lnv), 

+  hisv7 (lnv) , hisv8 (lnv) ,hisv9(lnv)  ,hisvlO(lnv) ,hisvll(lnv) , 

+  hisvl 2 (lnv) ,hisvl3 (lnv) ,hisvl4 (lnv) ,hisvl5 (lnv) ,hisvl6 (lnv) , 

+  hisvl 7 (lnv) , hisvl8 (lnv) ,hisvl9(lnv) ,hisv20(lnv) ,hisv21(lnv) , 

+  hisv22 (lnv) ,hisv23 (lnv) ,hisv24 (lnv) ,hisv25 (lnv) ,hisv26 (lnv) , 

+  hisv27 (lnv) , hisv2  8  (lnv) ,hisv29 (lnv) ,hisv30 (lnv) ,hisv31 (lnv) , 

+  hisv32 (lnv) ,hisv33 (lnv) ,hisv34 (lnv) ,hisv35 (lnv) ,hisv36 (lnv) , 

+  hisv37 (lnv) ,hisv38  (lnv)  ,hisv39 (lnv) ,hisv40 (lnv) ,hisv41 (lnv) , 

+  hisv42 (lnv) ,hisv43 (lnv) ,hisv44 (lnv) ,hisv45 (lnv) ,hisv46 (lnv) , 

+  hisv47 (lnv) ,hisv48  (lnv) ,hisv49 (lnv) ,hisv50 (lnv) ,hisv51 (lnv) , 

+  hisv52 (lnv) ,hisv53 (lnv) ,hisv54 (lnv) ,hisv55 (lnv) ,hisv56 (lnv) , 

+  hisv57 (lnv) ,hisv58 (lnv) , hisv59 (lnv) ,hisv60 (lnv) ,hisv61 (lnv) , 

+  hisv62 (lnv) fhisv63 (lnv) ,hisv64 (lnv) ,hisv65 (lnv) ,hisv66 (lnv) , 

+  hisv67 (lnv) ,hisv68  (lnv) ,hisv69 (lnv) ,hisv70 (lnv) ,hisv71 (lnv) , 

+  hisv72 (lnv) ,hisv73 (lnv) ,hisv74 (lnv) ,hisv75 (lnv) ,hisv76 (lnv) , 

+  hisv77 (lnv) ,hisv78 (lnv) ,hisv79 (lnv) ,hisv80 (lnv) ,hisv81 (lnv) , 

+  hisv82 (lnv) ,hisv83 (lnv) ,hisv84 (lnv) ,hisv85 (lnv) ,hisv86 (lnv) , 

+  hisv87 (lnv) ,hisv88 (lnv) ,hisv89 (lnv) ,hisv90 (lnv) ,hisv91 (lnv) , 

+  hisv92 (lnv) ,hisv93 (lnv) ,hisv94 (lnv) ,hisv95 (lnv) ,hisv96 (lnv) , 

+  hisv97 (lnv) ,hisv98 (lnv) ,hisv99 (lnv) 
common/ auxl 8 /dd (lnv) ,dfe(lnv) 

common/aux33/ixl (lnv) , ix2 (lnv)  ,  ix3 (lnv) , ix4 (lnv) , ix5 (lnv) , 

1  ix6 (lnv) , ix7 (lnv)  ,  ix8 (lnv) ,mxt (lnv) ,nmel 

common/ aux35/rhoa (lnv) ,cb(lnv) ,p(lnv) 
common/aux36/lf  t,  lit 
common/failu/sieu(lnv) , failu(lnv) 
common/sandb/isandb,nacbug, issflg 
common  / amain/ ipblank 
pointer ( ipblank , a  ( 1 ) ) 

dimension  cm(*) ,bqs (*) ,nhex(*) ,effs (lnv) 

REAL  N_E , N_S 
INTEGER  I,J 
INTEGER  NREG , NMAT 

DIMENSION  DREG (103 ,4,5) , DMAT (34 , 3 , 20)  , 

+  EE (3, 100)  ,  SIG_0E (3 , 100) ,N_E(3,100) , 

+  GG (3,100)  ,SIG_0S (3, 100) ,N_S(3,100) ,UEE(3,100) , 

+  ALPHA (6, 100) ,TH(100) ,ANGD(100) , DDT (100) , 
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+  HV (100+1) , ZV (100) , FR (5,3,100) , 

+  CTSIG{6, 6) , CTEPS (6,6) , CTSIGI (6 , 6) , CTEPSI (6 , 6) , 
+  EPSL (6) ,EPSG(6, 100) ,EPSP(6, 100) , 

+  ET (3 , 100) , GT (3,100) , 

+  CIJP (6, 6, 100) , CIJB (6, 6, 100) , ALPHAB (6, 100) , 

+  CCM  (6,6)  ,  CMNEW  (6,6)  ,  AB  (6,6)  , 

+  SIGL (6) ,RNU(3, 100) ,DEPSL(6) ,DSIGL(6) , 

+  TEMPI (6) ,TEMP2 (6) , TEMP3 (6) 

LOGICAL  KFIRST 
DATA  KFIRST/ . FALSE . / 

SAVE  KFIRST 

OPEN (UNIT=99 , FILE* 1 /home/ucn582/workarl/data . in ' 
+  , STATUS* 'OLD') 

C  OPEN (UNIT=1 , FILE* ' out ' , STATUS* ' old ' ) 

c 

c  write (*,*) 'cycles' ,  ncycle 

DO  45  1=1,6 
EPSL (I) =0.0 
TEMP2 (I) =0 . 0 
TEMP3 (I) =0 . 0 
DO  45  J=l, 6 

CMNEW  (I,  J)  =0.0 
CTSIG(i, j)=0.0 
CTEPS (i, j ) =0 .0 
CTSIGI (i, j ) =0 . 0 
CTEPSI (i, j ) =0 . 0 

45  CONTINUE 

DO  46  1=1,6 
DO  46  J=l, INPLY 
EPSG(I, J) =0 . 0 
EPSP(I, J)=0.0 

46  CONTINUE 
c 

c 

IF  (KFIRST)  GOTO  900 
c 

do  999  i=l,  6 
do  999  j-1,6 
CCM (i , j ) =0 . 0 
999  continue 

READ (99, *)  NREG 

READ (99, *)  NMAT 

write (1,*)  'NREG*' , NREG,  »  NMAT*', NMAT 

DO  12  IIA=1,NREG 

write (1,*)' - 

DO  22  IC=1 , 2 

READ (99, *)  (DREG(IC, IB, IIA) , IB=1,4) 
write(l,*)  (DREG (IC, IB, IIA)  ,IB=1,4) 

22  CONTINUE 

READ (99,*)  (DREG (3, IB, IIA) , IB=1, 4) 
write (1,*)  (DREG (IC, IB, IIA) ,IB=1,4) 

IND2=3 . 0+DREG (3 , 1, IIA) 

DO  32  IC-4.IND2 

READ (99, *)  (DREG(IC, IB, IIA) , IB=1,4) 
write (1,*)  (DREG (IC, IB, IIA) ,IB=1,4) 
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32  CONTINUE 

12  CONTINUE 

NLINES  =  34 

DO  10  IIA=1,NMAT 

write  (1,*)  1  MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM 1 
DO  10  I C=l, NLINES 

READ (99,*)  (DMAT ( IC, IB, IIA) ,16=1,3) 
write (1,*)  (DMAT (IC, IB, IIA) ,IB=1,3) 

10  CONTINUE 

C 

CLOSE (99) 

KFIRST  =  .TRUE. 

900  CONTINUE 

c 

mx=48* (mxt (1ft) -1) 

ym=cm(mx+l) 

pr=cm(mx+6) 

pass=cm(mx+16) 

fs=cm(mx+21) 

c  write (*,  *) ’ fs ' ,fs 

ss=cm(mx+2) 
g=ym/ (l.+pr) 
gdt=dtl*g 
gd2= . 5*gdt 

blk=-dtl*ym/ ((1.-2.  *pr) ) 
c 

do  910  i=lft,llt 
cb (i) =ss 

davg ( i ) =third*dd ( i ) 
p (i) =blk*davg (i) 

einc  (i)  =  (dl  (i)  *sigl  (i)  +d2  (i)  *sig2  (i)  +d3  (i)  *sig3  (i)  +d4  (i)  *sig4  (i) 
1  +d5 (i) *sig5 (i) +d 6 (i) *sig6 (i) +dd(i) *bqs (i) ) *dtl 

910  continue 
c 

c  write (*,*) 1 lft=  ' ,1ft, 'llt=  1 ,llt 

do  920  in=lft,llt 
c 
c 

IIA  =cm (mx+11) 

c  write (*,*) 'IIA1 , IIA 

FAIL_CRT  =  DREG  ( 1 , 2 , 1 1  A) 

NLOOP  =  DREG ( 1 , 3 , I I A) 

BETA  =  DREG (2,1, IIA) 

PHI  =  DREG  (2, 2,  IIA) 

SI  =  DREG  (2, 3,  IIA) 

NPLY  =  DREG ( 3 , 1 , I I A) 

DO  220  IC=1 , NPLY 

IDMAT  =  DREG ( IC+3 , 1 , IIA) 

c  write (*,*) 1 IDMAT= 1 , IDMAT 

DO  230  J=l, 3 

EE  ( J,  IC)  =  DMAT  (2,  J,  IDMAT) 

SIG_OE(J,IC)  =  DMAT  (3,  J,  IDMAT) 

N_E(J,IC)  =  DMAT  (4,  J,  IDMAT) 

RNU ( J, IC)  =  DMAT (5, J, IDMAT) 

GG  ( J,  IC)  =  DMAT  (6,  J,  IDMAT) 

SIG_OS(J,IC)  =  DMAT  (7  ,J,  IDMAT) 
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230 

241 


242 

220 

c 

C 


2000 

c 


c 

c 


11 


c 

c 


c 

123 

513 

C 


N_S(J,  IC) 
ALPHA (J, IC) 
UEE ( J, IC) 
CONTINUE 
DO  241  J=4 , 6 

ALPHA (J,IC) 
CONTINUE 

TH  (IC) 
ANGD(IC) 
DDT (IC) 
IROW 

DO  242  1=1,5 
DO  242  J=l, 3 

FR  (I ,  J,  IC) 
CONTINUE 
CONTINUE 


DMAT (8 , J, IDMAT) 
DMAT (9, J, IDMAT) 
DMAT (2, J, IDMAT) 


0.0 


DREG (IC+3 , 2 , IIA) 

DREG (IC+3 , 3 , IIA) 

DREG (IC+3, 4, IIA) 

10.0  +  (  ( FAIL_CRT- 1.0)  *  5.0  ) 


DMAT( (IROW-1) +I,J, IDMAT) 


write (*, *) 'hisv99 ' ,hisv99 (in) , 1  cycle' ,ncycle 
if  (hisv99 (in) .eq. 0 . 0)  goto  2000 
cnnl=hisv99 (in) *pass 
cnn=ncycle 

if  (cnn.lt.cnnl)  goto  1000 
hisv99 (in) =hisv99 (in) +1 . 0 

TEMPI (l)=hisvl (in) 

TEMPI (2) =hisv2 (in) 

TEMPI (3) =hisv3 (in) 

TEMPI (4) =hisv5 (in) 

TEMPI (5) =hisv6 (in) 

TEMPI (6) =hisv4 (in) 

CALL  KTRANSN  (BETA, PHI , SI,  CTSIG, CTEPS , CTSIGI, CTEPSI) 
CALL  KMATRIX  (CTEPS, TEMPI , EPSL, 6, 6 , 6, 6 , 6 , 6 , 1, 3 , IERR) 

DO  11  K=1,NPLY 

EPSG ( 1 , K) =EPSL (1) 

EPSG (2 , K) =EPSL (2) 

EPSG ( 3 , K) =EPSL (3 ) 

EPSG (4, K) =EPSL(4) 

EPSG (5, K) =EPSL (5) 

EPSG (6, K) =EPSL (6) 

CONTINUE 
BETA  =0.0 
PHI  =  0.0 
DO  123  K=1,NPLY 
SI=ANGD (K) 

write (1, *) '//////////  si////////', SI 

CALL  KTRANSN  (BETA, PHI , SI , CTSIG, CTEPS , CTSIGI , CTEPSI) 
CALL  KMATRIX  (CTEPS , EPSG (1, K) , EPSP (1, K) , 

+  6, 6, 6, 6, 6, 6, 1,3, IERR) 


CONTINUE 

format (6 (el4 . 8, lx, el4 . 8) ) 

CALL  KTANMOD  (NPLY, EPSP, EE, SIG_OE,N_E, GG, SIG_OS ,N_S, 
CALL  KCIJPRNLX  (NPLY, UEE,  ET,RNU, GT, CIJP) 


ET, GT) 
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c 

c  write(l,*) • *****************CUP************************* ' 

c  do  197  i=l, 6 

c  write (1,511)  (CIJP (i,  j , 1) , j=l, 6) 

cl97  continue 
c 

CALL  KBARS  (NPLY, ANGD, Cl JP, ALPHA, Cl JB , ALPHAB) 

CALL  KTHICK  (NPLY,TH,  HV.ZV.TOT) 

CALL  KCHOU  (NPLY, CIJB, TH, TOT,  CCM) 
c 

BETA  =  DREG (2,1, IIA) 

PHI  =  DREG (2, 2, IIA) 

SI  =  DREG  (2,3,1  IA) 

c  write (1,*)  'beta, phi, si ' .BETA, PHI, SI 

c  write(l,*) • *****************CUB************************* ’ 

c  do  194  1=1,6 

c  write  (1,511)  (CUB  (i,  j  ,  1)  ,  j=l,  6) 

cl94  continue 

c 

CALL  KTRANSN  (BETA, PHI , SI , CTSIG, CTEPS , CTSIGI , CTEPSI) 

CALL  KMATRIX  (CCM, CTEPS, AB, 6 , 6, 6, 6, 6 , 6, 6 , 3 , IERR) 
c 

c  write (1,*) ’ +++++++++++++++CTEPS++++++++++++++++ ' 

c  do  195  1=1,6 

c  write (1,511)  (CTEPS (i , j ) , j=l, 6) 

cl95  continue 

c 

CALL  KMATRIX  (CTSIGI , AB, CMNEW, 6, 6 , 6, 6 , 6 , 6 , 6, 3 , IERR) 
c 

hisv7 (in) =CMNEW (1,1) 
hisv8 (in) =CMNEW(1,2) 
hisv9 ( in) =CMNEW (1,3) 
hisvlO ( in) =CMNEW (1,4) 
hisvll (in) =CMNEW (1,5) 
hisvl2 (in) =CMNEW (1,6) 
hisvl3 (in) =CMNEW(2,2) 
hisvl4 (in) =CMNEW (2,3) 
hisvl5 ( in) =CMNEW (2,4) 
hisvl6 (in) =CMNEW (2,5) 
hisvl7 (in) =CMNEW (2,6) 
hisvl8 (in) =CMNEW (3,3) 
hisvl9 (in) =CMNEW (3,4) 
hisv20 (in) =CMNEW(3 , 5) 
hisv21 (in) =CMNEW (3 , 6) 
hisv22 (in) =CMNEW (4,4) 
hisv23 (in) =CMNEW (4,5) 
hisv24 ( in) =CMNEW (4,6) 
hisv25 ( in) =CMNEW (5,5) 
hisv26 (in) =CMNEW (5 , 6) 
hisv27 (in) =CMNEW (6,6) 
c 
c 

hisv28 (in) =CTEPS (1,  1) 
hisv29 (in) =CTEPS (1,2) 
hisv30 (in) =CTEPS (1,3) 
hisv31 (in) =CTEPS (1,4) 
hisv32 (in) =CTEPS (1,5) 
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hisv33 (in) =CTEPS (1,6) 
hisv34 (in) =CTEPS (2,2) 
hisv3 5 ( in) =CTEPS (2,3) 
hisv36 (in) =CTEPS (2,4) 
hisv37 ( in) =CTEPS (2,5) 
hisv38 (in) =CTEPS (2,6) 
hisv39 (in) =CTEPS (3,3) 
hisv40 (in) =CTEPS (3,4) 
hisv41 (in) =CTEPS (3,5) 
hisv42 (in) =CTEPS (3,6) 
hisv43 ( in) =CTEPS (4,4) 
hisv44 (in) =CTEPS (4,5) 
hisv45 (in) =CTEPS (4 , 6) 
hisv46 (in) =CTEPS (5,5) 
hisv47 (in) =CTEPS (5,6) 
hisv48 (in) =CTEPS (6,6) 
c 
c 

hisv4 9 ( in) =CTS IGI (1,1) 
hisvSO (in) =CTSIGI (1,2) 
hisv51 (in) =CTSIGI (1,3) 
hisv52 (in) =CTSIGI (1,4) 
hisv53 (in) =CTSIGI (1,5) 
hisv54 (in) =CTSIGI (1,6) 
hisv55 (in) =CTSIGI(2,2) 
hisv56 (in) =CTSIGI (2,3) 
hisv57 (in) =CTSIGI (2,4) 
hisv58 (in) =CTSIGI (2,5) 
hisv59 (in) =CTSIGI (2,6) 
hisv60 ( in) =CTSIGI (3,3) 
hisv61 ( in) =CTSIGI (3,4) 
hisv62 ( in) =CTSIGI (3,5) 
hisv63 (in) =CTSIGI (3,6) 
hisv64 (in) =CTSIGI (4,4) 
hisv65 (in) =CTSIGI (4,5) 
hisv66 (in) -CTSIGI (4,6) 
hisv67 (in) =CTSIGI (5,5) 
hisv68 (in) =CTSIGI (5,6) 
hisv69 (in) =CTSIGI (6,6) 
c 
c 

1000  continue 

c 

c 

CMNEW (1, 1) =hisv7 (in) 
CMNEW(1,2) =hisv8 (in) 
CMNEW (1 , 3) =hisv9 (in) 
CMNEW (1,4) =hisvl0 (in) 
CMNEW (1, 5) =hisvll (in) 
CMNEW (1,6) =hisvl2 (in) 
CMNEW (2,2) =hisvl3 (in) 
CMNEW (2,3) =hisvl4 (in) 
CMNEW (2,4) =hisvl5 (in) 
CMNEW (2,5) =hisvl6 (in) 
CMNEW (2 , 6) =hisvl7 (in) 
CMNEW (3,3) =hisvl8 (in) 
CMNEW ( 3 , 4 ) =hi s vl 9 ( in) 
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CMNEW (3,5) =hisv20 (in) 
CMNEW (3,6) =hisv21 ( in) 
CMNEW (4,4) =hisv22 (in) 
CMNEW (4,5) =hisv23 (in) 
CMNEW (4,6) =hisv24 ( in) 
CMNEW (5, 5) =hisv25 (in) 
CMNEW  (5, 6)  =hisv26  (in) 
CMNEW  (6, 6)  =hisv27  (in) 
c 

CMNEW  (2,1)  =hisv8  ( in) 
CMNEW  (3,1)  =hisv9  ( in) 
CMNEW  (3,2)  =hisvl4  (in) 
CMNEW (4,1) =hisvlO (in) 
CMNEW (4,2) =hi s vl5 ( in) 
CMNEW (4,3) =hisvl9 (in) 
CMNEW (5,1) =hisvll ( in) 
CMNEW (5,2) =hisvl6 (in) 
CMNEW (5,3) =hisv20 (in) 
CMNEW (5,4) =hisv23 (in) 
CMNEW (6,1) =hisvl2 (in) 
CMNEW (6,2) =hisvl7 (in) 
CMNEW  (6,3)  =hi  s  v2 1  ( in) 
CMNEW (6,4) =hisv24 (in) 
CMNEW (6, 5) =hisv26 (in) 
c 
c 

CTEPS (1,1) =hisv28 (in) 
CTEPS (1,2) =hi sv2 9 ( in) 
CTEPS (1, 3) =hisv30 (in) 
CTEPS (1,4) =hi sv3 1 ( in) 
CTEPS (1, 5) =hisv32 (in) 
CTEPS (1,6) =hisv33 (in) 
CTEPS (2,2) =hi sv3 4 ( in) 
CTEPS (2,3) =hi sv3 5 ( in) 
CTEPS (2,4) =hisv36 (in) 
CTEPS (2,5) =hisv37 (in) 
CTEPS (2,6) =hi sv3 8 ( in) 
CTEPS (3,3) =hisv39 (in) 
CTEPS (3,4) =hisv40 (in) 
CTEPS (3,5) =hisv41 (in) 
CTEPS (3,6) =hisv42 (in) 
CTEPS (4,4) =hi sv4 3 ( in) 
CTEPS (4,5) =hi sv44 (in) 
CTEPS (4,6) =hisv45 (in) 
CTEPS (5,5) =hi sv4  6 ( in) 
CTEPS (5,6) =hisv47 (in) 
CTEPS (6,6) =hi sv4  8 ( in) 
c 

CTEPS (2,1) =hisv29 (in) 
CTEPS ( 3 , 1 ) =hi s v3  0 ( in) 
CTEPS (3,2) =hisv35 (in) 
CTEPS (4,1) =hisv31 (in) 
CTEPS (4,2) =hisv36 (in) 
CTEPS (4 , 3) =hisv40 (in) 
CTEPS (5,1) =hisv32 (in) 
CTEPS (5,2) =hi sv3 7 ( in) 
CTEPS (5,3) =hisv41 (in) 
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CTEPS (5,4) =hisv44 ( in) 

CTEPS (6,1) =hisv33 (in) 

CTEPS (6,2) =hisv3 8 ( in) 

CTEPS (6,3) =hisv42 (in) 

CTEPS (6,4) =hisv45 (in) 

CTEPS (6,5) =hisv47 (in) 
c 
c 

CTSIGI (1, 1) =hisv49 (in) 

CTSIGI (1,2)  =hisv50  (in) 

CTSIGI (1,3) =hisv51 (in) 

CTSIGI (1,4) =hisv52 (in) 

CTSIGI (1,5) =hisv53 (in) 

CTSIGI (1,6) =hi sv54 (in) 

CTSIGI (2 , 2 ) =hisv55 ( in) 

CTSIGI (2,3) =hisv56 (in) 

CTSIGI (2,4) =hisv57 (in) 

CTSIGI (2,5) =hisv58 (in) 

CTSIGI (2,6) =hisv59 (in) 

CTSIGI (3,3) =hisv60 (in) 

CTSIGI (3,4) =hisv61 (in) 

CTSIGI (3,5) =hisv62 (in) 

CTSIGI (3,6) =hisv63 (in) 

CTSIGI (4,4) =hisv64 (in) 

CTSIGI (4,5) =hisv65 (in) 

CTSIGI (4,6) =hisv66 (in) 

CTS IGI (5,5) =hi s v6  7  ( in) 

CTSIGI (5,6) =hisv68 (in) 

CTSIGI (6,6) =hisv69 (in) 
c 

CTSIGI (2,1) =hi sv50 (in) 

CTSIGI (3,1) =hisv51 (in) 

CTSIGI (3,2) =hisv56 (in) 

CTSIGI (4 , 1)  =hi sv52 (in) 

CTSIGI (4,2) =hisv57  ( in) 

CTSIGI (4,3) =hisv61 (in) 

CTSIGI (5, l)=hisv53 (in) 

CTSIGI (5,2) =hisv58 (in) 

CTSIGI (5,3) =hisv62 (in) 

CTSIGI (5,4) =hisv65 (in) 

CTSIGI (6, l)=hisv54 (in) 

CTSIGI (6,2) =hisv59 (in) 

CTSIGI (6,3) =hisv63 (in) 

CTSIGI (6,4) =hisv66 (in) 

CTSIGI (6,5)  =hisv68 (in) 
c 

TEMP2 (l)=dl (in) *dtl 
TEMP2 (2) =d2 (in) *dtl 
TEMP2 (3)=d3 (in) *dtl 
TEMP2 (4) =d4 (in) *dtl 
TEMP2 (5) =d5 (in) *dtl 
TEMP2 (6) =d6 (in) *dtl 
c 

CALL  KMATRIX  (CTEPS, TEMP2 ,DEPSL, 6, 6, 6, 6, 6 , 6, 1, 3  ,  IERR) 
CALL  KMATRIX  (CMNEW, DEPSL, DSIGL, 6 , 6 , 6 , 6 , 6 , 6, 1, 3 , IERR) 
CALL  KMATRIX  (CTSIGI, DSIGL, TEMP3, 6, 6, 6, 6, 6, 6, 1, 3,  IERR) 
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sigl (in) =sigl (in) +TEMP3 (1) 
sig2 (in) =sig2 (in) +TEMP3 (2) 
sig3 (in) =sig3 (in) +TEMP3 (3) 
sig4 (in) =sig4 (in) +TEMP3 (4) 
sig5 (in) =sig5 (in) +TEMP3 (5) 
sig6 (in) =sig6 (in) +TEMP3 (6) 
c 

hisvl (in) =hisvl (in) +temp2 (1) 
hisv2 (in) =hisv2 (in) +temp2 (2) 
hisv3 (in)=hisv3 (in)+temp2 (3) 
hisv4 (in) =hisv4 (in) +temp2 (6) 
hisv5 (in) =hisv5 (in) +temp2 (4) 
hisv6 (in) =hisv6 (in) +temp2 (5) 
c 
c 

c  write (1,*)  *  strain:  stress:1 

c  write (1,*)  hisvl (in) , sigl (in) 

c  write (1,*)  hisv2 (in) , sig2 (in) 

c  write (1,*)  hisv3 (in) , sig3 (in) 

c  write (1,*)  hisv4 (in) , sig4 (in) 

c  write (1,*)  hisv5 (in) , sig5 (in) 

c  write (1,*)  hisv6 (in) , sig6 (in) 

c  write(l,*)  'Stiffness  Matrix;  ddsdde:1 

do  193  i=l , 6 

c  write (1,511)  (CMNEW(i, j) , j-1,6) 

193  continue 

c  write  (1,  *)  '&&&&&:&&&&&&&&&&&:&&&  end  &&&&&&&&&&&&&&&&&&&&& 1 

511  format (6 (el4 . 8, lx) ) 

c 

c 

c 

c 

920  continue 
c 

c  f s=0 .20 

do  80  i=lft,llt 

ielmtc=lnv* (ndum-1) +i 
ielmtc=nhex  (ielmtc) 

ef fs (i) =sqrt (2 . * (hisvl (i) **2+hisv2 (i) **2+hisv3 (i) **2 
1  +0.5* (hisv4 (i) **2+hisv5 (i) **2+hisv6 (i) **2) ) /3 . ) 

c  write ( * , * ) i , nhex ( i ) , ielmtc,eff s (i) 

if  (failu(i) .eq.0.0)  go  to  80 
if (isandb.eq.l.and.effs (i) .gt.fs)  then 
write  (*  ,5040)  ncycle, ielmtc 
failu(i)=0.0 
sigl (i) =0 . 
sig2 (i) =0 . 
sig3(i)=0. 
sig4 (i) =0 . 
sig5 (i) =0 . 
sig6 (i) =0 . 
endif 

cine (i) — (dl (x) *oigl (i) id2 (i) *oig2 (i) id3 (i) +oig3 (i) idd (i) *cig4 (i> 

1  +d5 (i) *sig5 (i) +d6 (i) *sig6 (i) ) *dtl+einc (i) 

80  continue 
c 

5040  format (/'  AT  THIS  CYCLE  ' ,i8,'  THE  FOLLOWING  ELEMENT  FAILED  ',i8/) 
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